library(sasld)

# -------------------------------------------------------------------
# ATTENZIONE: la funzione sample e' stata modificata 
# a partire da R 3.6.0 per cui i "vecchi risultati"
# non posso essere riprodotti di default
# http://hostingwin.unitn.it/micciolo/PeM/01-rng.pdf
# -------------------------------------------------------------------
# se si vogliono riprodurre i risultati dei Laboratori
# con le versioni pił recenti di R va eseguita la seguente istruzione
# -------------------------------------------------------------------
RNGkind(sample.kind="Rounding")
# -------------------------------------------------------------------


# 8.1 ---------------------------------------------------------------
set.seed(246135)
out <- simIC.prop(p=0.5,n=50,nrep=100)
summary(out)
summary(out, conf.level=0.9)
summary(out, conf.level=0.99)
plot(out)
plot(out, conf.level=0.9)

set.seed(246135)
out <- simIC.prop(p=0.5,n=50,nrep=100000)
summary(out, conf.level=0.99)

set.seed(246135)
out <- simIC.prop(p=0.5,n=50,nrep=100)
plot(out, conf.level=0.9)
tmp <- summary(out)
showIntervals(tmp)
# -------------------------------------------------------------------




# 8.2 ---------------------------------------------------------------
set.seed(162534)
out <- simIC.mean(mean=170,sd=10,n=20,nrep=100)
summary(out)
plot(out)
plot(out, conf.level=0.9)
#
tmp <- summary(out)
showIntervals(tmp)

# pochissime osservazioni
set.seed(162534)
out <- simIC.mean(mean=170,sd=10,n=4,nrep=10000)
summary(out)


# uniforme
set.seed(162534)
out <- simIC.mean(n=5,nrep=10000,family="unif",min=10,max=20)
summary(out)

set.seed(162534)
out <- simIC.mean(n=10,nrep=10000,family="unif",min=10,max=20)
summary(out)


# esponenziale
set.seed(123456)
out <- simIC.mean(n=10,nrep=10000,family="exp",rate=2)
summary(out)
hist(out)
tmp <- out$out[,1]; qqnorm(tmp); qqline(tmp)

set.seed(123456)
out <- simIC.mean(n=100,nrep=10000,family="exp",rate=2)
summary(out)
hist(out)
tmp <- out$out[,1]; qqnorm(tmp); qqline(tmp)

set.seed(123456)
out <- simIC.mean(n=1000,nrep=10000,family="exp",rate=2)
summary(out)
hist(out)
tmp <- out$out[,1]; qqnorm(tmp); qqline(tmp)


# dati "reali" (distribuzione empirica)
set.seed(123456)
out <- simIC.mean(n=10,nrep=10000,family="real",data=c(0:9))
summary(out)
hist(out)
tmp <- out$out[,1]; qqnorm(tmp); qqline(tmp)
# -------------------------------------------------------------------




# 8.3 ---------------------------------------------------------------
set.seed(123456)
out <- simIC.prop(p=0.1,n=10,nrep=100)
summary(out)

set.seed(123456)
out <- simIC.prop(p=0.1,n=10,nrep=100000)
summary(out)

set.seed(123456)
out <- simIC.prop(p=0.1,n=10,nrep=100)
plot(out)


#
set.seed(654321)
out <- simIC.prop(p=0.1,n=150,nrep=100000)
summary(out)
hist(out)
boxplot(out)
tmp <- out$out[,1]; qqnorm(tmp); qqline(tmp)

set.seed(654321)
out <- simIC.prop(p=0.1,n=500,nrep=10000)
summary(out)
hist(out)
tmp <- out$out[,1]; qqnorm(tmp); qqline(tmp)

set.seed(654321)
out <- simIC.prop(p=0.1,n=5000,nrep=10000)
summary(out)
hist(out)
tmp <- out$out[,1]; qqnorm(tmp); qqline(tmp)
# -------------------------------------------------------------------




# 8.4 ---------------------------------------------------------------
ph <- 18/20
es <- (ph*(1-ph))/20; es <-sqrt(es); es
me <- 1.96*es; me
c(ph-me,ph+me)

dbinom(18,20,0.8)+dbinom(19,20,0.8)+dbinom(20,20,0.8)
dbinom(18,20,0.7)+dbinom(19,20,0.7)+dbinom(20,20,0.7)
dbinom(18,20,0.65)+dbinom(19,20,0.65)+dbinom(20,20,0.65)

p <- 0.65; dbinom(18,20,p)+dbinom(19,20,p)+dbinom(20,20,p)

p <- seq(0.65,0.7,0.001)
y <- dbinom(18,20,p)+dbinom(19,20,p)+dbinom(20,20,p)
plot(p,y)
abline(h=0.025,lty=3)
tmp <- cbind(p,y)
p <- 0.683; dbinom(18,20,p)+dbinom(19,20,p)+dbinom(20,20,p)

# estremo superiore
sum(dbinom(c(0:18),20,0.95))
sum(dbinom(c(0:18),20,0.99))
p <- seq(0.97,0.999,0.001)
y <- numeric(length(p))
for (i in 1:length(p)) y[i] <- sum(dbinom(c(0:18),20,p[i]))
plot(p,y)
abline(h=0.025,lty=3)
tmp <- cbind(p,y)
p <- 0.9877; sum(dbinom(c(0:18),20,p))

#
binom.test(18,20)
p <- 0.6830173; dbinom(18,20,p)+dbinom(19,20,p)+dbinom(20,20,p)
p <- 0.9876515; sum(dbinom(c(0:18),20,p))
# -------------------------------------------------------------------




# 8.5 ---------------------------------------------------------------
w <- c(160.2,160.8,161.4,162.0,160.8,162.0,162.0,161.8,161.6,161.8)
c(mean(w),sd(w))
sort(w)

set.seed(654321)
sample(w,10,replace=TRUE)

set.seed(654321)
sort(sample(w,10,replace=TRUE))
sd(sample(w,10,replace=TRUE))
sd(sample(w,10,replace=TRUE))


#
set.seed(123456)
x <- sample(w,1000000,replace=TRUE)
B <- matrix(x,nrow=100000)
y <- apply(B,1,sd)
hist(y,xlab="Deviazione standard campionaria",ylab="Frequenza",main="")
box()
quantile(y,prob=c(0.025,0.975),type=1)
tmp <- sort(y); c(tmp[2500],tmp[97500])
quantile(y,prob=c(0.05,0.95),type=1)


# media (letture di peso)
y <- apply(B,1,mean)
hist(y,xlab="Media campionaria",ylab="Frequenza",main="")
box()
quantile(y,prob=c(0.025,0.975),type=1)
#
n <- length(w)
mn <- mean(w)
ds <- sd(w)
es <- ds/sqrt(n)
tcr <- qt(0.025,(n-1),lower.tail=FALSE)
me <- tcr*es
c(mn-me,mn+me)
#
tcr <- qt(0.05,(n-1),lower.tail=FALSE)
me <- tcr*es
c(mn-me,mn+me)
quantile(y,prob=c(0.05,0.95),type=1)
# -------------------------------------------------------------------